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Abstract 


We address the problem of planning collision-free paths 
for multiple agents using optimization methods known 
as proximal algorithms. Recently this approach was ex¬ 
plored in [Bento et al.| j2013^ , which demonstrated its 
ease of parallelization and decentralization, the speed 
with which the algorithms generate good quality solu¬ 
tions, and its ability to incorporate different proximal 
operators, each ensuring that paths satisfy a desired 
property. Unfortunately, the operators derived only ap¬ 
ply to paths in 2D and require that any intermediate 
waypoints we might want agents to follow be preas¬ 
signed to specific agents, limiting their range of applica¬ 
bility. In this paper we resolve these limitations. We in¬ 
troduce new operators to deal with agents moving in ar¬ 
bitrary dimensions that are faster to compute than their 
2D predecessors and we introduce landmarks, space- 
time positions that are automatically assigned to the set 
of agents under different optimality criteria. Finally, we 
report the performance of the new operators in several 
numerical experiment^ 

1 Introduction 


In this paper we provide a novel set of algorithmic build¬ 
ing blocks (proximal operators) to plan paths for a system 
of multiple independent robots that need to move optimally 
across a set of locations and avoid collisions with obstacles 
and each other. This problem is crucial in applications in¬ 
volving automated storage, exploration and surveillance. 

Even if each robot has few degrees of freedom, the joint 


system is complex and this problem is hard to solve (Reif 
|1979[popcroft, Schwartz, and Sharir 1984) . We can divide 
existing algorithms for this problem into global planners, 
if they find collision-free beginning-to-end paths connecting 
two desired conhgurations, or local planners, if they find 
short collision-free paths that move the system only a bit 
closer to the hnal configuration. 

We briefly review two of the most rigorous approaches. 


Random sampling methods, first introduced in ( [Kavraki and 
[Latombe 1994[[Kavraki et al. 1996| l, are applicable to global 
planning and explore the space of possible robot configura¬ 
tions with discrete structures. The rapidly-exploring random 
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tree algorithm (RRT; La Valle and Kuffner 200 l| l, is guar¬ 
anteed to asymptotically find feasible solutions with high- 


probabilty while the RRT* algorithm (Karaman and Fraz- 


zoli 2010 1 asymptotically finds the optimal solution. How¬ 


ever, their convergence rate degrades as the dimension of 
the joint configuration space increases, as when considering 
multiple robots, and they cannot easily find solutions where 
robots move in tight spaces. In addition, even approximately 
solving some simple problems requires many samples, e.g., 
approximating a shortest path solution for a single robot re¬ 
quired to move between two points with no obstacles (|Kara- 


man and Frazzoli 2010 see Fig. 1). These methods explore 


a continuous space using discrete structures and are differ¬ 
ent from methods that only consider agents that move on a 
graph with no concern about their volume or dynamics, e.g. 
(Standley and Korf [Sharon et al. 2013) . 

An optimization-based approach has been used by sev¬ 
eral authors, including IMellinger, Kushleyev, and Kumar 


(20121, who formulate global planning as a mixed-integer 
quadratic problem and, for up to four robots, solve it using 
branch and bound techniques. Sequential convex program¬ 
ming was used in (lAugugliaro, Schoellig, and D’Andrea 


|2012j l to efficiently obtain good local optima for global plan¬ 
ning up to twelve robots. State-of-the-art optimization-based 
algorithms for local planning typically have real-time per¬ 
formance and are based on the velocity-obstacle (VO) idea 
of ( Fiorini and Shiller 1998| ), which greedily plans paths 
only a few seconds into the future and then re-plans. These 
methods scale to hundreds of robots. Unlike sampling al¬ 
gorithms, optimization-based methods easily find solutions 
where robots move tightly together and solve simple prob¬ 
lems very fast. However, they do not perform as well in 
problems involving robots in complex mazes. 

Our work builds on the work of jBento et al.| ( |2013] l, which 
formulates multi-agent path planning as a large non-convex 
optimization problem and uses proximal algorithms to solve 
it. More specifically, the authors use the Alternating Direc¬ 
tion Method of Multipliers (ADMM) and the variant intro¬ 
duced in Derbinsky et al. (2013 1 called the Three Weight 
Algorithm (TWA). These are iterative algorithms that do not 
access the objective function directly but indirectly through 
multiple (simple) algorithmic blocks called proximal opera¬ 
tors, one per function-term in the objective. At each iteration 
these operators can be applied independently and so both the 














































TWA and the ADMM are easily parallelized. A brief expla¬ 
nation of this optimization formulation is given in Section 
1^ A self-contained explanation about proximal algorithms 
is found in |Parikh and Boyd (|20131 l and a good review on 
the ADMM is |Boyd et al.| ( |2()1 l| l. In general the ADMM 
and the TWA are not guaranteed to converge for non-convex 
problems. There is some work on solving non-convex prob¬ 
lems using proximal algorithms with guarantees (see|Udell| 
and Boyd 2014 and references in|Parikh and Boyd|2013|l but 


the settings considered are not applicable to the optimization 
problem at hand. Nonetheless, the empirical results in|Bento| 


et al. (2013 1 are very satisfactory. For global planning, their 


algorithm scales to many more robots than other optimiza¬ 
tion based methods and finds better solutions than VO-based 
methods. Their method also can be implemented in the (use¬ 
ful) form of a decentralized message-passing scheme and 
new proximal operators can be easily added or removed to 
account for different aspects of planning, such as, richer dy¬ 
namic and obstacle constraints. 

The main contributions of Bento et al. ( 2013| l are the prox¬ 
imal operators that enforce no robot-robot collisions and no 
robot-wall collisions. These operators involve solving a non¬ 
trivial problem with an infinite number of constraints and a 
finite number of variables, also known as a semi-infinite pro¬ 
gramming problem (SIP). The authors solve this SIP prob¬ 
lem only for robots moving in 2D by means of a mechanical 
analogy, which unfortunately excludes applications in 3D 
such as those involving fleets of unmanned aerial vehicles 
(UAVs) or autonomous underwater vehicles (AUVs). An¬ 
other limitation of their work is that it does not allow robots 
to automatically select waypoint positions from a set of ref¬ 
erence positions. This is required, for example, in problems 
involving robots in formations (Bahceci, Soysal, and Sabin] 
20031. In|Bento et al. (20131, any reference position must be 


pre-assigned to a specific robot. 

In this paper we propose a solution to these limitations. 
Our contributions are (i) we rigorously prove that the SIP 
problem involved in collision proximal operators can be re¬ 
duced to solving a single-constraint non-convex problem 
that we solve explicitly in arbitrary dimensions and numeri¬ 
cally show our novel approach is substantially faster for 2D 
than Bento et ar| (2013 1 and (ii) we derive new proximal op¬ 
erators that automatically assign agents to a subset of refer¬ 
ence positions and penalize non-optimal assignments. Our 
contributions have an impact beyond path planning prob¬ 
lems. Other applications in robotics, computer vision or 
CAD that can be tackled via large optimization problems 
involving collision constraints or the optimal assignment of 


objects to positions (e.g. Kuffner et al. 2002 


Kass|198^ Andreev, Pavisic, and Raspopovic||2001[) might 


Witkin and 


benefit from our new building blocks (cf. Section]^. 

While there is an extensive literature on how to solve SIP 
problems (see |Stein ( |2012 1 for a good review), as far as we 
know, previous methods are either too general and, when ap¬ 
plied to our problem, computationally more expensive than 
our approach, or too restrictive and thus not applicable. 

Finally, we clarify that our paper is not so much about 


solved critical limitations. However, our numerical results 
and supplementary video do confirm that the framework pro¬ 
duces very good results, although there are no guarantees 
that the method avoids local minima. 

2 Background 


Here we review the formulation of [Bento et al.| P013] l of 
path planning as an optimization problem, explain what 
proximal operators are, and explain their connection to solv¬ 
ing this optimization problem. 

We have p spherical agents in of radius Our 

objective is to find collision-free paths {a:i(f)}iG[p],tG[o,T] 
for all agents between their specified initial positions 
at time 0 and specified final positions 
at time T. In the simplest case, we divide time in intervals 
of equal length and the path {xi{t)}J^Q of agent i G [p] is 
parametrized by a set of break-points {a;i(s)}g^g such that 
Xi{t) =Xi{s) for t = sT/rj and all s. Between break-points 
agents have zero-acceleration. We discuss the practical im¬ 
pact of this assumption in Appendix |A] 

We express global path planning as an optimization prob¬ 
lem with an objective function that is a large sum of simple 
cost functions. Each function accounts f or a different aspec t 
of the problem. Using similar notation tojBento et al.|(|2013|l, 
we need to minimize the objective function 

i i 

fi°j (Xi (s) , Xi (s-fl) , Xj (s) , Xj (s-l-1)) -fE fr'{xi{s) , x^ (s-Fl)) 

i,s 

+ E fw"iXi{s),Xi{s-i-l)) +^ff'’^{Xi{s),Xi{s-\-l),Xi{s-i-2)). 

W,i,S i,S 

The function /‘^°** prevents agent-agent collisions: it is 
zero if ||aa:i(s)-|-(l —a)a:i(s-|-l) —(aa;j(s)-|-(l —a)a:j(s-|- 
l))l! > U + ^j for all a G [0,1] and infinity otherwise. The 
ywaii function prevents agents from colliding with obstacles: 
it is zero if \\axi{s) -f (1 — Q;)a:i(s -f 1) — y|| > for all 
a G [0,1], y G W, where W is a set of points defining an 
obstacle, and is infinity otherwise. In |Bento et al.| ( 20T3] l, W 
is a line between two points x l and in the plane and the 
summation across a set of obstacles. The functions 

Z''®' and 7*^“ impose restrictions on the velocities and direc¬ 
tion changes of paths. The function /p°* imposes boundary 
conditions: it is zero if Xi (0) = (or if Xi (rj) = x^™’) and 

infinity otherwise. The authors also re-implement the local 
path planning method of [Alonso-Mora et al.|p013| l, based 
on velocity obstacles ( Fiorini and Shiller 1998| ), by solving 
an optimization algorithm similar to IF 

[Bento et ar] ( |2013| l solve (0 using the TWA, a variation of 
the ADMM. The ADMM is an iterative algorithm that mini¬ 
mizes objectives that are a sum of many different functions. 
The ADMM is guaranteed to solve convex problems, but, 
empirically, the ADMM (and the TWA) can find good fea¬ 


showing the merits of the framework used in Bento et al. 
(2013[ a point already made), as it is about overcoming un¬ 


sible solutions for large non-convex problems (Derbinsky et 
[al. 2013[|Bento et al. 20131 . 

Loosely speaking, the ADMM proceeds by passing mes¬ 
sages back and forth between two types of blocks: proximal 
operators and consensus operators. First, each function in 









































































the objective is queried separately by its associated prox¬ 
imal operator to estimate the optimal value of the variables 
the function depends on. For example, the proximal operator 
associated with produces estimates for optimal value of 
xi(s), a:: 2 (s), a:i(s -f 1) and X 2 (s + 1). These estimates are 
then sent to the consensus operators. Second, a consensus 
value for each variable is produced by its associated con¬ 
sensus operator by combining all the received different es¬ 
timates for the values of the variable that the proximal op¬ 
erators produced. For example, the proximal operators as¬ 
sociated with fi °2 and give two different estimates for 
the optimal value of a;i(s) and the consensus operator asso¬ 
ciated with a;i(s) needs to combine them into a single es¬ 
timate. The consensus estimates produced by the consensus 
operators are then communicated back to and used by the 
proximal operators to produce new estimates, and the cycle 
is repeated until convergence. See Appendix [B] for an illus¬ 
tration of the blocks that solve a problem for two agents. 

It is important to be more specific here. Consider a func¬ 
tion f{x) in the objective. From the consensus value for its 
variables x, the corresponding consensus nodes form con¬ 
sensus messages n that are sent to the proximal operator as¬ 
sociated with /. The proximal operator then estimates the 
optimal value for cc as a tradeoff between a solution that is 
close to the minimizer of / and one that is close to the con¬ 
sensus information in n ( [Parikh and Boyd 201 3| l, 

a: G argmin/(a;') +-||a;'- n|p, ( 2 ) 

x' Z 


where we use G instead of = to indicate that, for a non- 
convex function /, the operator might be one-to-many, in 
which case some extra tie-breaking rule needs to be imple¬ 
mented. The variable p is a free parameter of the ADMM 
that controls this tradeoff and whose value affects its perfor¬ 
mance. In the TWA the performance is improved by dynam¬ 
ically assigning to the p’s of the different proximal operators 
values in {0, const., oo} (cf. Appendix [Q. 

We emphasize that the implementation of these proximal 
operators is the crucial inner-loop step of the ADMM/TWA. 
For example, when / = /™* takes a quadratic (kinetic en¬ 
ergy) form, the operator 0 has a simple closed-form ex¬ 
pression. However, for / = or f — the operator 
involves solving a SIP problem. In Section|^we explain how 
to compute these operators more efficiently and in a more 
general setting than in Bento et ar](|2013|l. 


3 No-collision proximal operator 


Here we study the proximal operator associated with the 
function that ensures there is no collision between two 
agents of radius r and r' that move between two consecutive 
break-points. We distinguish the variables associated to the 
two agents using ’ and distinguish the variables associated 
to the two break-points using _ and “, respectively. For con¬ 
creteness, just imagine, for example, that x = Xi{0),x' = 
2 ^ 2 ( 0 ), ai = a;i(l) and x' = a; 2 (l) and think of n, n', n' and 
n' as the associated received consensus messages. Following 


the operator associated to /' 


coll 


P, 

r,T' ,x,x' Z 


f + ^l 

" ^ 2 ' 


f + ^ 

II -r 2 


outputs the minimizer of 

2 


\x —n 


f+ -I 

II T 2 I 


X —n 


s.t. \\a{x — x') + {l — a){x — x')\\ > r-|-r', for all aG[0,1]. (3) 


Our most important contribution here is an efficient pro¬ 
cedure to solve the above semi-infinite programming prob¬ 
lem for agents in arbitrary dimensions ^ reducing it to a 
max-min problem. Concretely, Theorem [^below shows that 

t is essentially equivalent to the ‘most costly’ of the prob- 
as in the following family of single-constraint problems 
parametrized by a, 


■ Pu 
mm =1 

x.x' ,x,x' Z 


l|2 , P||_ —I 

II +2ll^-«l 


I^ + -Il 

I -r 2 II 


l|2 , P 11—/ _/ 

II +2 11*-^ 


S.t. \\a{x — x') + {1 — a){x — x')\\ > r + r'. 


(4) 


Since problem Q has a simple closed-form solution, we can 
solve (|^ faster than in Bento et al. (20131 for 2D objects. We 
support this claim with numerical results in Section]^ In the 
supplementary video we use our new operator to do planning 
in 3D and, for illustration purposes, in 4D. 

Theorem 1. /f|ja(n —n')-I-(1 — a)(n —n')|j 7 ^ 0, f/zen 
has a unique minimizer, x* (a), and if this condition holds 
for a = a* G arg max^/gjo,!] h{a'), where 2hf(a) is the 
minimum value of 0 . then x*{a*) is also a minimizer of 
([^. In addition, (/'||Q!(ri — n') -f (1 — a)(n — n') || 7 ^ 0, then 


h{a) = max < 0, 


(r -I- r') — ||aAn -|- (1 — a)Ar 


(1 - ay/p 

and the unique minimizer of 0 is 

X* = n — 'yp{a^ An + a{l — a)An), 
x'* = n + ')p {a^An -\- a(l — a) An), 

X* =n — 7 p((l — a)aAn -|- (1 — a)^An), 

X* =n-\- 7 p^((l — a)aAn -|- (1 — a)^An), 

where p = {p~^ + p = {p~^ + p'~^) 

_ 2A \ _ 

T — T I o \ P,-.,2 / I ^.\2 / > A — 


(5) 


( 6 ) 

(7) 

( 8 ) 

(9) 

-C-l 


l+ 2 A(aVP+(l-a)Vp) 

An = n — n! and An = n — n'. 


Ho:) 

f(r+r') /p+(l-ay / f 


Remark 2. Under a few conditions, we can use Theorem 
\I\tofind one solution to problem (|^ by solving the simpler 
problem (|^ for a special value of a. In numerical imple¬ 
mentations however, the conditions of Theorem^are easy to 
satisfy, and the x*{a*) obtained is the unique minimizer of 
problem (|^. We sketch why this is the case in Appendix^\ 

In a nutshell, to find one solution to we simply find a* 
by maximizing 0 and then minimize 0 using ( 0-0 with 
a = a*. We can carry both steps efficiently, as shown in 
Section|5] The intuition behind Theorem[T]is that if we solve 
the optimization problem 0 for the ‘worst’ constraint (the 
a* that gives largest minimum value), then the solution also 
satisfies all other constraints, that is, it holds for all other 
a G [0,1]. We make this precise in the following general 
lemma that we use to prove Theorem We denote by di the 
derivative of a function with respect to the variable. The 
proof of this Lemma is in Appendix]^ and that of Theorem 
[T]is in Appendix]^ 

Lemma 3. Let A be a convex set in K, p : x ^ 

K, (x, a) 1 —^ g{x, a), be convex in a and continuously dif¬ 
ferentiable in (x, a) and let f : —>■ M, x 1 —>■ /(x). 

















be continuously differentiable and have a unique minimizer. 
For every a & A, let h{a) denote the minimum value of 
min 3 ,.g( 2 . q,)>q f{x) and if the minimum is attained by some 
feasible point let this be denoted by x*{a). Under these 
conditions, if a* G arginaxQg_4 and if x* (a) exists 
around a neighborhood of a* in A and is differentiable at 
a*, and if dig(x*(a*), a*) f 0, then a;*(Q;*) minimizes 

^i^a::3(3;,Q;)>OVae:,4 /{^)' 


Other collision operators 

Using similar ideas to those just described, we now explain 
how to efficiently extend to higher dimensions the wall- 
agent collision operator that [Bento et al. (2013 i introduced. 
In the supplementary video we use these operators for path 
planning with obstacles in 3D. 

To avoid a collision between agent 1, of radius r, and a 
line between points yi,y 2 & we include the following 
constraint in the overall optimization problem; ||axi{s) -I- 
(l-a)a;i(s + l)-(/3j/i + (l-/ 3 )?/ 2 )|| > r for all a, /3 S 
[0,1] and all s + 1 € [ry + 1]. This constraint is associated 
with the proximal operator that receives (n,n') and finds 
{x,x) that minimizes |||x — n||^ + |||x — n|p subject to 
||ax+(l-a)x-(/3yi + (l-/ 3 ) 2 / 2 )|| >r,foralla,^ G [0,1]. 
Using ideas very similar to those behind Theorem and 
Lemma we solve this problem for dimensions strictly 
greater than two by maximizing over a,P £ [ 0 , 1 ] the min¬ 
imum value the single-constraint version of the problem. In 
fact, it is easy to generalize Lemma to ^ C and the 
single-constraint version of this optimization problem can be 
obtained from (ffl by replacing nf and W with /3r/i-|-(l—/ 3 )j/ 2 , 
and letting p',p —>■ oo. Thus, we use © and ([^-([^ under 
this replacement and limit to generalize the line-agent colli¬ 
sion proximal operator of jBento et al.| ( |20T3l ) to dimensions 
greater than two. We can also use the same operator to avoid 
collisions between agents and a line of thickness u, by re¬ 
placing r with u + r. 

Unfortunately, we cannot implement a proximal operator 
to avoid collisions between an agent and the convex enve¬ 
lope of an arbitrary set of points yi,y 2 , ■■■,yq by maximiz¬ 
ing over a, /3i,..., Pq-i G [0,1] the minimum of the single¬ 
constraint problem obtained from (0 after replacing rf and 
n' with l3iyi+...+j3q-2yq-2+{y-lii----l3q-i)yq, and letting 
p',~p' oo. We can only do so when d > q, otherwise we 
observe that the condition dig{x*{a*),a*) f 0 of Lemma|^ 
does not hold and x* (a*) is not feasible for the original SIP 
problem. In particular, we cannot directly apply our max- 
min approach to re-derive the line-agent collision operator 
for agents in 2D but only for dimensions > 3. When d < q, 
we believe that a similar but more complicated principle can 
be applied to solve the original SIP problem. Our intuition 
from a few examples is that this involves considering differ¬ 
ent portions of the space A separately, computing extremal 
points instead of maximizing and minimizing and choosing 
the best feasible solution among these. We will explore this 
further in future work. 


Speeding up computations 

The computational bottleneck for our collision operators is 
maximizing Here we describe two scenarios, denoted 


as trivial and easy, when we avoid this expensive step to 
improve performance. 

First notice that one can readily check whether x = n 
is a trivial feasible solution. If it is yes, it must be optimal, 
because it has 0 cost, and the operator can return it as the 
optimal solution. This is the case if the segment from An = 
n — n' to An = n — W does not intersect the sphere of 
radius r-l-r' centered at zero, which is equivalent to ||aAn-|- 
(1 — ajAnjj > r + r' with a = max{l, min{ 0 , a'}} and 
a' = Aff {An — An)/|| An — Anp. 

The second easy case is a shortcut to directly determine if 
the maximizer of max„g[o h{a) is either 0 or 1. We start 
by noting that empirically h has at most one extreme point 
in [0,1] (the curious reader can convince him/herself of this 
by plotting h{a) for different values of An and An). This 
being the case, if i9i/i(0) > 0 and dih{l) > 0 then a* = 1 
and if dih{0) < 0 and dih{l) < 0 then a* = 0. Evalu¬ 
ating two derivatives of h is much easier than maximizing 
h and can save computation time. In particular, dih{0) — 
C{—{r -\- r') + ||An|| + (An^(An — An)/||An||)) and 
dih{l) = C'{{r + r') — || Anjj -f (An'l'(An — An)/|| Anjj)) 
for constants C, C > 0. 

If these cases do not hold, we cannot avoid maximizing 
a scenario we denote as expensive. In Sectionj^we pro¬ 
file how often each scenario occurs in practice and the cor¬ 
responding gain in speed. 


Local path plauuiug 

The optimization problem 0 finds beginning-to-end 
collision-free paths for all agents simultaneously. This is 
called global path planning. It is also possible to solve path 
planning greedily by solving a sequence of small optimiza¬ 
tion problems, i.e. local path planning. Each of these prob¬ 
lems plans the path of all agents for the next t seconds such 
that, as a group, they get as close as possible to their fi¬ 
nal desired positions. This is done, for example, in Eiorini 
and Shiller ( 1998|l and followup work ( jAlonso-Mora et ^ 
2012a[ jAlonso-Mora et al. 2013] l. The authors in [Bento ^ 
al.[ ( [2013| l solve these small optimization problems using a 
special case of the no-collision operator we study in Section 
and show this approach is computationally competitive 
with the results in |Alonso-Mora et al.|P013) . Therefore, our 
results also extend this line of research on local path plan¬ 
ning to arbitrary dimensions and improve solving-times even 
further. See Section]^ for details on these improvements in 
speed. 

4 Landmark proximal operator 

In this section we introduce the concept of landmarks that, 
automatically and jointly, (i) produce reference points in 
space-time that, as a group, agents should try to visit, (ii) 
produce a good assignment between these reference points 
and the agents, and (iii) produce collision-free paths for the 
agents that are trying to visit points assigned to them. 

Points (i) to (iii) are essential, for example, to formation 
control in multi-robot systems and autonomous surveillance 
or search (Bahceci, Soysal, and Sabin 20031, and are also 


related to the problem of assigning tasks to robots, if the 
tasks are seen as groups of points to visit (Michael et al. 



























2008 1 . Many works focus on only one of these points or treat 
them in isolation. One application where points (i) to (iii) 
are considered, although separately, is the problem of using 
color-changing robots as pixels in a display (|Alonso-Mora 


et al. 2012bl |Alonso-Mora et al. 2012ct |Ratti and Frazzoli' 


2010 1 . The pixel-robots arrangement is planned frame-by- 
frame and does not automatically guarantee that the same 
image part is represented by the same robots across frames, 
creating visual artifacts. Our landmark formalism allows us 
to penalize these situations. 

We introduce landmarks as extra terms in the objective 
function we now explain how to compute their associ¬ 
ated proximal operators. Consider a set of landmark trajecto¬ 
ries and, to each trajectory j, assign 

a cost Cj > 0 , which is the cost of ignoring the entire land¬ 
mark trajectory. In addition, to each landmark yj{s) G 
that is assigned to an agent, assign a penalty Cj (s) > 0 for 
deviating from yAs). Landmark trajectories extend the ob¬ 
jective function dTjl by adding to it the following term 


s = Sinit 


where the variable tjj indicates which agent should follow 
trajectory j. If aj = *, that means trajectory j is unassigned. 
Each trajectory can be assigned to at most one agent and 
vice-versa, which it must follow throughout its duration. So 
we have aj G [p] U {*} as well as the condition that if 
aj = aj! then either j = f or aj = *. We optimize the 
overall objective function over x and a. Note that it is not 
equally important to follow every point in the trajectory. For 
example, by setting some c’s equal to zero we can effectively 
deal with trajectories of different lengths, different begin¬ 
nings and ends, and even trajectories with holes. By setting 
some of the c’s equal to infinity we impose that, if the trajec¬ 
tory is followed, it must be followed exactly. In ( flOl i we use 
the Euclidean metric but other distances can be considered, 
even non-convex ones, as long as the resulting proximal op¬ 
erators are easy to compute. Finally, notice that, a priori, we 
do not need {t/j(s)} to describe collision free trajectories. 
The other terms in the overall objective function will try to 
enforce no-collision constraints and additional dynamic con¬ 
straints. Of course, if we try to satisfy an unreasonable set of 
path specifications, the ADMM or TWA might not converge. 

The proximal operator associated to term receives as 
input {ni(s)} and outputs {a;*(s)} where i G [p], Sinit < 
s < Send and {a:*(s)} minimizes 


®end 

min X! 0(s)ll®-T, (s) Cj 

s = Sinji j-(Tj=* 

P 5 end 

+ 5 Z X! - ni{s)f. ( 11 ) 

i=l s = Sjjjj{ 


The variables cr’s are used only internally in the computa¬ 
tion of the proximal operator because they are not shared 
with other terms in the overall objective function. The above 
proximal operator can be efficiently computed as follows. 
We first optimize over the x’s as a function of a and 
then we optimize the resulting expression over the cr’s. If 


we optimize over the cc’s we obtain J^j where, if 
aj = *, ujj^if = Cj and, if aj = i A ujjj = 


min.E:=s. 

■'^end pi Cj ( 


s=Si„i, Cj (s j 


Xiis) - Pjf -f f ||xi(s) - ni(s)|p = 
2 c^{s)+^ ll^t(s) - . The last equality fol¬ 

lows from solving a simple quadratic problem. We can op¬ 
timize over the cr’s by solving a linear assignment prob¬ 
lem with cost matrix w, which can be done, for example, 
using Hungarian method of |Kuhn| (|1955[), using more ad- 


vanced methods such as those after Goldberg and Tarjan 
(|1988 1 , or using scalable but sub-optimal algorithms as in 
[Bertseka^ ( |1988 1 . Once an optimal cr* is found, the output 
of the operator can be computed as follows. If i is such that 
{j : cr*=i} = 0 then x* (s) = n^{s) for all Sinit < s < Send 
and if i is such that i = aj for some j G [m] then x* (s) = 
{pini{s)+2Ci{s)yj{s))/{2Ci{s)+Pt) for all Sink < S < Send- 

The term ( [T 0 | ) corresponds to a set of trajectories between 
break-points s = Sink and s = Send for which the different 
agents must compete, that is, each agent can follow at most 
one trajectory. We might however want to allow an agent 
to be assigned to and cover multiple landmark trajectories. 
One immediate way of doing so is by adding more terms 
of the form ( [TOl i to the overall objective function such that 
the term has all its trajectories within the inter- 

['^inin'^end]’ ^^'d different intervals for different k’s are 
disjoint. However, just doing this does not allow us to im¬ 
pose a constraint like the following; “the trajectory in 
the set corresponding to the interval [sLj, Sg^^J must be cov¬ 
ered by the same agent as the the {j'Y^ trajectory in the 

set corresponding to the interval ■^end^^^]” "f® do so 

we need to impose the additional constraint that some of 
the variables across different terms of the form ( [T 0 | ) 
are the same, e.g. in the previous example, 

Since the variables a’s can now be shared across different 
terms, the proximal operator o needs to change. Now it 
receives as input a set of values {ni(s)}s,i and and 

outputs a set of values s and {tr*}j that minimize 

In the expression above, {aj}j and {n'are both vectors 
of length p -I- 1, where the last component encodes for no 
assignment and aj must be binary with only one 1 entry. 
For example, if p = 5 and a 2 = [0, 0,1,0,0,0] we mean 
that the second trajectory is assigned to the third agent, or 
if (74 = [0,0,0, 0, 0, 1 ] we mean that the fourth trajectory is 
not assigned to any agent. However, n' can have real values 
and several nonzero components. 

We also solve the problem above by first optimizing over 
X and then over a. Optimizing over x we obtain J2j ’ 


where = Wj,*-f ^ || [ 0 , ... 0 , 1 , 0 ,..., 0 ] - n' 

4ll<f+lll-n 


Given the cost matrix w, we find the optimal cr* by solving 
a linear assignment problem. Given a*, we compute the op¬ 
timal X* using exactly the same expressions as for 0- 


— Wi.i+ 


tIIS 


J 

/ ||2 




+l- 2 n 


/(i) 


(k) 

Finally, to include constraints of the kind a) = a 


1^') 

























we add to the objective a term that takes the value inhnity 
whenever the constraint is violated and zero otherwise. This 
term is associated with a proximal operator that receives as 
input n' = and n'-, = and 

outputs (o-*,cr*,) G argmin^^.^^^., ^\\(Jj-n'j\\'^ + ^\\aj'- 
n' / p. Again aj and aj' are binary vectors of length p + 1 
with exactly one non-zero entry. The solution has the form 
a* = a*, = [0, 0,..., 0,1, 0, ...0] where the 1 is in position 

i = arg maxig[pj pjn> + p^,n>, . 

5 Numerical experiments 


cores decrease time (we see ~ 2x with 8 cores). In Figure 
[T]-(c) we show that the paths found when the TWA solves 
CONFl in 3D over 1000 random initializations are not very 
different and seem to be good (in terms of objective value). 


(a) 

CONF1-2D 


CONFl-3D*%^ § 


In CONFl, antipodal 
agents around a circle or 
sphere exchange position 



Number of agents, p 


(c) 



We gathered all results with a Java implementation of the 
ADMM and the TWA as described in [Bento et al.| ( |2013t 
see Appendix jC]! using JDK7 and Ubuntu vl2.04 run on a 
desktop machine with 2.4GHz cores. 

We hrst compare the speed of the implementation of the 
collision operator as described in this paper, which we shall 
refer to as “NEW,” with the implementation described in 
[Bento et al.j ( [2013[ ), which we denote “OLD.” We run the 
TWA using OLD on the 2D scenario called “CONLl” in 
[Bento et'ar] ( [2013] l with p = 8 agents of radius r = 0.918, 
equally spaced around a circle of radius i? = 3, each re¬ 
quired to exchange position with the corresponding antipo¬ 
dal agent (cf. Lig. [T]-(a)). While running the TWA using 
OLD, we record the trace of all n variables input into the 
OLD operators. We compare the execution speed of OLD 
and NEW on this trace of inputs, after segmenting the n 
variables into trivial, easy, or expensive according to Lor 
global planning, the distribution of trivial, easy, expensive 
inputs is {0.814,0.001,0.185}. Although the expensive in¬ 
puts are infrequent, the total wall-clock time that NEW takes 
to process them is 76 msec compared to 54 msec to pro¬ 
cess all trivial and easy inputs. By comparison, OLD takes 
a total time of 551 msec on the expensive inputs and so our 
new implementation yields an average speedup of 7.25 x on 
the inputs that are most expensive to process. Similarly, we 
collect the trace of the n variables input into the collision 
operator when using the local planning method described 
in [Bento et al.[ ( [2013[ l on this same scenario. We observe 
a distribution of the trivial, easy, expensive inputs equal to 
(0.597,0.121,0.282}, we get a total time spent in the easy 
and trivial cases of 340 msec for NEW and a total time spent 
in the expensive cases of 2802 msec for NEW and 24157 
msec for OLD. This is an average speedup of 8.62x on the 
expensive inputs. Lor other scenarios, we observe similar 
speedup on the expensive inputs, although scenarios easier 
than CONLl normally have fewer expensive inputs. E.g., if 
the initial and final positions are chosen at random instead of 
according to CONLl, this distribution is (0.968, 0, 0.032}. 

Ligure [T]-(b) shows the convergence time for instances of 
CONLl in 3D (see Lig. [^(a)) using NEW for a different 
number of agents using both the ADMM and the TWA. We 
recall that OLD cannot be applied to agents in 3D. Our re¬ 
sults are similar to those in [Bento et al.[ ( |^13[ l for 2D: (i) 
convergence time seems to grow polynomially with p; (ii) 
the TWA is faster than the ADMM; and, (iii) the proxi¬ 
mal operators lend themselves to parallelism, and thus added 


Ligure 1: (a) CONL1-2D & 3D; (b) Convergence time for 
CONL1-3D for a varying number of cores and agents; (c) 
Empirical distribution of the objective over 1000 random ini¬ 
tializations of TWA for CONL1-3D. 

In the supplementary movie we demonstrate the use of the 
landmark operators. Lirst we show the use of these operators 
on six toy problems involving two agents and four landmark 
trajectories where we can use intuition to determine if the so¬ 
lutions found are good or bad. We solve these six scenarios 
using the ADMM with 100 different random initializations 
to avoid local minima and reliably find very good solutions. 
With 1 core it always takes less than 3 seconds to converge 
and typically less than 1 second. We also solve a more com¬ 
plex problem involving 10 agents and about 100 landmarks 
whose solution is a ‘movie’ where the different robots act 
as pixels. With our landmark operators we do not have to 
pre-assign the robots to the pixels in each frame. 

6 Conclusion 

We introduced two novel proximal operators that allow 
the use of proximal algorithms to plan paths for agents in 
3D, 4D, etc. and also to automatically assign waypoints to 
agents. The growing interest in coordinating large swarms 
of quadcopters in formation, for example, illustrates the im¬ 
portance of both extensions. Lor agents in 2D, our collision 
operator is substantially faster than its predecessor. In par¬ 
ticular, it leads to an implementation of the velocity-obstacle 
local planning method that is faster than its implementation 
in both [Alonso-Mora et al.| ( |2013| l and [Bento et al.[ ( |2013| l. 
The impact of our work goes beyond path planning. We are 
currently working on two other projects that use our results. 
One is related to visual tracking of multiple non-colliding 
large objects and the other is related to the optimal design of 
layouts, such as for electronic circuits. In the first, the speed 
of the new no-collision operator is crucial to achieve real¬ 
time performance and in the second we apply Lemma to 
derive no-collision operators for non-circular objects. 

The proximal algorithms used can get stuck in local min¬ 
ima, although empirically we find good solutions even for 
hard instances with very few or no random re-initializations. 
Luture work might explore improving robustness, possibly 
by adding a simple method to start the TWA or the ADMM 
from a ‘good’ initial point. Linally, it would be valuable to 
implement wall-agent collision proximal operators that are 
more general than what we describe in Section]^ perhaps 
by exploring other methods to solve SIP problems. 
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Appendix for “Proximal Operators for 
Multi-Agent Path Planning” 


A A comment on the impact of the 
assumption of piece-wise linear paths in 
practice 


A direct application of the approach of [Bento et al. (2013 1 
can result in very large accelerations at the break-points. It 
is however not hard to overcome this apparent limitation in 
practical applications. We now explain one way of doing it. 
First notice that by increasing the number of break-points 
we can obtain trajectories arbitrarily close to smooth trajec¬ 
tories with finite acceleration everywhere. Since in practice 
it is not efficient to work with a very large number of break¬ 
points, we can keep the number of break-points at a rea¬ 
sonable level and increase the effective radius of the robots. 
This would allow us to fit a polynomial through the break¬ 
points and obtain smooth trajectories that are never distant 
from the piece-wise linear paths by more than the difference 
between the true robots radii and their effective radii. Us¬ 
ing this approach, we would obtain a set of non-colliding 
finite-acceleration trajectories. We can also impose specific 
maximum-acceleration constrains if, at the same time, we 
restrict the maximum permitted change of velocity at break¬ 
points using additional proximal operators. 


B An illustration of the two kinds of blocks 
used by the ADMM/TWA 

As explained in Section]^ the ADMM/TWA is an iterative 
scheme that alternates between (i) producing different esti¬ 
mates of the optimal value of the variables each function in 
the objective depends on and (ii) producing consensus val¬ 
ues from the different estimates that pertain the same vari¬ 
able. The blocks that produce the estimates we call proxi¬ 
mal operators and the blocks that produce consensus values 
we call consensus operators. The proximal operator blocks 
only receive messages from the consensus blocks (and send 
estimates back to them) and the consensus blocks only re¬ 
ceive estimates from the proximal operators (and send con¬ 
sensus messages back to the proximal operators). Hence, 
the ADMM iteration scheme can be interpreted as messages 
passing back and forth along the edges of a bipartite graph. 

We now illustrate this. Imagine that we want to com¬ 
pute non-colliding paths for two agents and that trajecto¬ 
ries are parametrized by three break-points. In this case 
the optimization problem Q has six variables, namely, 
a;i(0),a:i(l),a;i(2) for agent 1 and a; 2 ( 0 ),a; 2 (l),a; 2 ( 2 ) for 
agent 2. See Figure[^(top). 

The ADMM has one consensus operator associated with 
the position of each agent at each break-point (blue blocks 
with “=’ sign on it) and has one proximal operator associated 
with each function in the objective (red blocks with function 
names on it). In our small example, we have two no-collision 
operators. One ensures there are no-collisions between the 
segment connecting a;i(0) and a;i(l) and the segment con¬ 
necting 2 : 2 ( 0 ) and 2 : 2 ( 1 ). The other acts similarly on a:i(l), 
a:i(2), 2 : 2 ( 1 ) and 2 : 2 ( 2 ). We also have four velocity oper¬ 
ators that penalize trajectories in which the path segments 








Figure 2: (Top) Variables in the problem; (Center) Graph 
of connections between proximal operators and consensus 
nodes; (Bottom) Proximal operators’ input and output val¬ 
ues. 


have large velocities. Finally, we have four position proxi¬ 
mal operators that enforce agent 1 and 2 to start and finish 
their paths at specified locations. See Figure [^(center). 

The crucial part of implementing the ADMM or the TWA 
is the construction of the proximal operators. The proximal 
operators receive consensus messages n from the consensus 
nodes and produce estimates for the values of the variables 
of the function associated to them. For example, at each iter¬ 
ation, the left-most no-collision proximal operator in Figure 
l^(center) receives messages n from the consensus nodes as¬ 
sociated to the variables 2 : 1 ( 0 ), 2 : 1 ( 1 ), 2 : 2 ( 0 ) and 2 : 2 ( 1 ) and 
produces new estimates for their optimal value. The center- 


































top consensus operator receives fours estimates for a;i(l), 
from two no-collision proximal operators and two velocity 
proximal operators, and produces a single estimate its opti¬ 
mal value. See Figure [^(bottom). 


C A comment on our implementation of the 

TWA 

Implementing the TWA requires computing proximal oper- 
ato rs and specifying, at every iteration, what Derbinsky et 
^(2013 I calls the outgoing weights, , of each proximal 
operator. In the TWA there are also incoming weights, ^, 
which correspond to the p’s that appear in the definition of 
all our proximal operators, but their update scheme at every 
iteration is fixed ( Derbinsky et al. 2013[ Section 4.1). 

For the collision operator and landmark operator we intro¬ 
duce in this paper, we compute the outgoing weights using 
the same principle as in ( [Bento et al. 2013 l. Namely, if an 
operator maps a set of input variables {ni,n2, ■■■,nk) to a 
set of output variables {xi,X 2 , ■■■,Xk) then, if rii ^ Xi we 
set = 1 otherwise, when Ui = Xi, we set pi = 0. In 
other words, if a variable is unchanged by the operators, the 
outgoing weight associated with it should be zero, otherwise 
it is 1 . 


D More details about Remark m 

The next three points sketch why Remark[^is true. 

First, if the p’s are all positive, the objective function of 
problem 0 is strictly convex and we can add the constraint 
\\{x,xf,x,x')\\ < M, for M large enough, without chang¬ 
ing its minimum value. The new extended set of constraints 
amounts to the intersection of closed sets with a compact set 
and by the continuity of the objective function and the ex¬ 
treme value theorem it follows that the problem has a mini- 
mizer. 

Second, problem ([^ can be interpreted as resolving the 
following conflict. “Agent 1 wants to move from n at time 
0 to n at time 1 and agent 2 wants to move from n' at time 
0 to n' at time 1 , however, if both move in a straight-line, 
they collide. How can we minimally perturb their initial and 
final reference positions so they they avoid collision?” If the 
vectors (n, 0), (n',0), (n, l),(n',l) lie in the same three- 
dimensional plane there can be ambiguity on how to min¬ 
imally perturb the agents’ initial and final positions; agent 
I’s reference positions can either move ‘up’ and agent 2 ’s 
‘down’ or agent I’s reference positions ‘down’ and agent 
2 ’s ‘up’ (‘up’ and ‘down’ relative to the plane defined by 
the vectors (n, 0), (n',0), (n, l),(n', 1)). In numerical im¬ 
plementations however, it almost never happens that (n, 0), 
(n', 0), (n, l),(n', 1 ) lie in the same plane and, in fact, this 
can be avoided by adding a very small amount of random 
noise to the n’s before solving problem ([^. 

Third, one can show from the continuity the objective 
function of problem 0 , the continuous-differentiability of 
g{x,a) = ||Q;(a; — x') -f (1 — a){x — — (r -|- r'Y 

and the fact that the level sets of p(x, a) as a function of x 
never have ‘flat’ sections that h{a) is a continuous function. 
Since [0,1] is compact, it follows that there always exists 
an a*. Also, for the purpose of a numerical implementation. 


we can consider ||a*(n — n') -|- (1 — a*){n — n')|j ^ 0. In 
fact, this can be avoided by adding a very small amount of 
random noise to the n’s before solving problem (|^. There¬ 
fore, for practical purposes, we can consider that for each 
a* there exists a unique minimizer to problem (|^. Finally, 
a careful inspection of 0 shows that if there exists a such 
that h(a) > 0 then a* is unique and if h{a) = 0 for all a 
then 0-(0 always give a; = n. In short, for the purpose of a 
numerical implementation, we always And a unique x*{a*) 
and because, in practice, as argued in the two points above, 
problem 0 can be considered to have only one solution, it 
follows that this unique x*{a*) is the unique minimizer of 
problem 0 . 


E Proof of Lemma |3] 

We need to consider two separate cases. 

In the first case we assume that x*{a*) is such 
that g(x*{a*),a*) > 0. This implies that a;*(Q;*) 

is a minimizer of min^; f{x), which implies that 
^^^x:g(x,a -)>0 fix) = miuj,/(x), which im¬ 
plies that min,j,.g(a;,Q,)>o/(a:) > min^fix) = 

fix) = maxa/g^ min,r:3(a:.a')>o/(a:) > 
min 2 ,,g(,r,a)>o fix)- In other words, min,r:g(a:.a)>o fix) = 
fix*ia*)) for all a € -4. Since fix) has a unique mini¬ 
mizer we have that x* (a*) must be feasible for the problem 
^^^x:g(x,a)>o fix), which implies that gix*ia*),a) > 0 
for all a G A. In other words, x* (a*) is feasible point of 
nnmx:g(^x,a)>oyaGA fix) and attains the smallest possible 
objective value, hence it minimizes it. 

In the second case we assume that gix*ia*),a*) = 0. To 
finish the proof it suffices to show that vd2gix*ia*), a*) > 
0 for all u € K such that a* + v G .4. To see this, we 
first notice that, if this is the case, then, for all a G A, 
we have by convexity that ( 7 (x*(a*), a) > gix*ia*),a*) + 
ia — a*)d2gix*ia*),a*) = (a — a*)d2gix*ia*), a*) > 
0 , which implies that x*ia*) is a feasible point for 
the problem nimx-g(x,a)>o,\faGA fix)- Secondly, we no¬ 
tice that /(x*(a*)) > min^.g(^_„)>o,vaG^/(a;) > 

maxaeA'^^^x-.g{x,a)>o fix) = fix*ia*))- This implies 
that x*ia*) minimizes min 3 ,.g( 3 ._a)>o,vaGyt fix)- 

We now show that vd 2 gix*ia*), a*) > 0 for all u S K 
such that a* -\- V G A- 

First we notice that since / is differentiable and x*(a) 
exists and is differentiable at x*(a*), then, by the chain 
rule, /i(a) is differentiable at a*- In addition, we notice 
that if a* maximizes h then, if u S K is such that a* -\- 

V G A, the directional derivative of h in the direction of 

V evaluated at a* must be non-positive. In other words, 
vdihia*) = vdifix*ia*))^dix*ia*) < 0. Second, we 
notice that in a small neighborhood around a*, x*(a) ex¬ 
ists and is continuous (because x*(a*) is differentiable) 
which, by the continuous-differentiability of g and the fact 
that digix*ia*),a*) 7 ^ 0 , implies that digix*ia),a) 7 ^ 0 
in this neighborhood. Therefore, in a small neighborhood 
around a*, the problem min^-.gi^^, (j)>o /(x) has a single in¬ 
equality constraint and digix*ia),a) 7 ^ 0 which implies 
that x*(a), which we are assuming exists in this neighbor¬ 
hood, is a feasible regular point and satisfies the first-order 











necessary optimality conditions (Bertsekas 1999f 


dif{x*{a)) + Xdig{x*{a),a) = 0, 
Xg{x*{a),a) = 0, 
g{x*{a),a) > 0 , 
A < 0. 


( 12 ) 

(13) 

(14) 

(15) 


Now, we take the directional derivative of ( [T4l l with respect 
to a in the direction v evaluated at (a;*(a*), a*) and obtain 

vdig{x*(a*), a*ydix* (a*) + vd 2 g{x*{a*),a*) > 0 . 

(16) 

At the same time, by computing the inner product 
between (EH and dix* evaluated at a* and mul¬ 
tiplying by V we obtain vdif{x*{a*)ydix*{a*) + 
Xvdig{x*{a*),a*ydix*(a*) = 0. But we have al¬ 

ready proved that vdif(x*{a*))^dix*{a*) < 0 therefore 
Xvdig{x*{a*),a*)^dix*{a*) > 0. Now recall that we are 
assuming g(x*{a*),a*) = 0 therefore, A < 0. It thus fol¬ 
lows that vdig{x*{a*),a*)''^dix*{a*) < 0 and from ( [T^ 
we conclude that vd 2 g{x*{a*), a*) > 0 . 


F Proof of Theorem [T] 

We first observe that the solutions to problem © and 0 
remain the same if we replace |ia(x — x') -f (1 — a)(x — 
^011 ^ ?’-l-r’'by \\a{x — xf) + {\ — a){x — x')\\‘^ — {r+r'Y > 
0. We prove the theorem with this replacement in mind. 

We first prove the theorem assuming that we have proved 
that that expressions 0 to 0 hold when ||Q!(zr — g!) -f (1 — 
a){n — n')|| ^ 0. This is then proved last. 

To prove the theorem we consider two separate cases. 

In the first case we consider r + r' < minQ,g[o^i] ||a(n — 
n') -f (1 — a){n — n')\\. This implies that the minimum of 
problem 0 is 0 and, because the p’s are positive, there is a 
unique minimizer which equals (n, n', n, n'). In this case we 
also have h{a) = 0 for all a, which implies that the solution 
of problem 0 is x*{a) = n for all a. Therefore, for any 
optimal a* for which \\a*{n — ri') + {1 — a*){n — n')\\ 7 ^ 0 , 
we have that the minimizer of problem 0 is equal to the 
unique solution of problem 0 . Hence the theorem is true in 
this case. 

Now we assume that r + r' > min^gjo^i] ||cr(n — n') -f 
(1 — a){ri — n')|j. This implies that there exists an a for 
which the right-hand-side of 0 is positive and for which 
\\a{n — n') -f (1 — a){n — ^|j ^ 0. This implies that, 
there exists an a for which hia) > 0. Therefore, for any 
optimal a* it must the case that h{a*) > 0. If ||cr*(n — 
n') -f (1 — a*){n — n')\\ 7 ^ 0 then 0 holds around a 
neighborhood of a* and in this neighborhood h{a) > 0 . 
Therefore, in this neighborhood x*{a) exists, is unique, and 
is differentiable. Now we notice that the objective func¬ 
tion of problem 0 is continuously-differentiable, and that 
g{x, a) = ||a(x — x') + (1 — a)(x — x')\\‘^ — {r + r')^ is 
continuously-differentiable in {x, a) and convex in a. If it is 
true that c)ip(a;*(a*), a*) 7 ^ 0 then we can apply Lemmaj^ 
with A = [ 0 , 1 ] and conclude that x*{a*) is also a solution 

'Also known as also known as Kamsh-Kuhn-Tucker (KKT) 
conditions. 


of problem 0. Hence the theorem will be true in this case 
as well. 

We now show that in this case, indeed, 
dig{x*{a*),a*) 7 ^ 0. First we notice that dig = 

{dxg,dxg,dx'g,dx>g) = 2 ||a*(x* - x*') -I- (1 - a*){x* - 
lE*^)||(a*,l — a*, —a*, —1 + a*). We now recall that, as 
explain above, if ||a*(ri — n') -f (1 — a*){n — n')\\ 7 ^ 0 
then h{a*) > 0. Therefore, we conclude that 

r + r' = ||a*(x* — x*') -f (1 — a*){x* — because 
otherwise — x*') + (1 — Q!*)(x* — x*')|| < r + r' 

implies that the constraint of problem 0 for a = a* is 
inactive which implies that the solution must he x = n with 
objective value 0 which contradicts the fact that h{a*) > 0 . 
Hence, dig = (r -f r'){a*, 1 — a*, —a*, —1 + a*), which 
is always non-zero, and proves that the theorem is true in 
this case as well. 

Finally, we now prove that that expressions 0 to 0 hold 
when \\a{n — n') -f (1 — a)(n — fi')|| 7 ^ 0. This amounts 
to a relatively long calculus computation and is written in 
Appendix [G| 


G Computation of minimum value and 
minimizer of problem (|0 

We do not solve problem 0 but instead solve the equivalent 
(more smooth) problem 

min_Jk-nf-f (17) 

x,x' ,x,x Z Z 



s.t. ||q!(x — a;') + (1 — a){x — x')|p — (r -f r'Y > 0. 

To begin, we notice that if 

\\a{n — n') + (1 — a){n — n')|| > (r -f r') 

then the constraint is inactive and the minimizer of ( fT7] i is 
{n,n',n,n') with minimum value 0. At the same time, if 
\\a{n — it/) + {1 — a){n — n')|| > (r + r') we have that 
h{a) = 0 and the minimizer obtained from 0-0 is also 
(n, n', n, n'). Therefore, we only need to show that equa¬ 
tions 0 and 0-0 hold in the case when the constraint is 
active, which corresponds to the case when \\a{ri — n') + 
(1 — a)(n — n')\\ < (r -f r'). 

To do so, we first introduce a few block variables (writ¬ 
ten in boldface) and express the above problem in a shorter 
form. Namely, we define x = (x, x',3;, x') S n = 

(n, n', n, n') S and a = {a, —a, (1 —a), —(1 —a)) S 
and D = diag(p, p', p, p') G and, rewrite ( [T7| ) as 

min -tr{(x — n)^I?(x — n)} 

X 2 (18) 

s.t. ||a^x||^ — (r -f > 0 . 

Then we notice that it is necessary that the solutions to this 
problem are among the points that satisfy the KKT condi¬ 
tions. Namely, those points that satisfy 


D{x — n) -f 2otv = 0 

(19) 

V = A(q:^x) 

( 20 ) 

lkll/|A = r + r' 

( 21 ) 




where A 7 ^ 0 is the Lagrange multiplier associated to the 
problem’s constraint and is non-zero because we are assum¬ 
ing the constraint is active. In the rest of the proof we show 
that there are only two points that satisfy the KKT condi¬ 
tions and show that, between them, the one that corresponds 
to the global optimum satisfies (|^ and ([^-( 1 ^. 

We first write the two equations even more compactly as 



We claim that, if 1 
block matrix 


2Xa:^D-^ a ^ 0, the inverse of the 


\D 


a. 


a 


t _ 


1/A 


(23) 


2{D-^ - 


2D~^Xaa^ D 


tn-i 


l-|-2AQ;tD- 
2Xa^D-^ 
l+2XaW~^cy. 


2XD-^a 

l-|-2AatD-ic 

-A 

T+2XoJTF^ 


■ (24) 


To prove this we could use the formula for the inverse of 
a block matrix. Instead, and much more simply, we simply 
compute the product of ( |24l i and ( |2^ and show it equals the 
identify. It is immediate to see thaUhe block diagonal entries 
of the resulting product are indeed identity matrices. Since 
both matrices are symmetric, all that is left to check is that 
one of the non-diagonal block entries is zero. Indeed, 


(at) 2 O-^- 


1 -H 2 AatD-i( 


+ (-l/A) 


2 AatD- 


^1-1- 2AatD-ia, 

_ 2a^D-\l + 2AatD-ia) - ia''D-^Xcxa^D-^ - 2a^D-^ 
~ 1 -h 2AatD-ia 

4Xa^ D~^ oc^ D~^ a — 4a^ D~^ Xaa^ D~^ 

^ l + 2XaW-^ct 

_ (at7:)-ia)(4AatD-i-4AatT>-i) _ 

“ l + 2AatD-ia ^ ' 

We now solve the linear system ( |2^ by multiplying both 
sides by the inverse matrix ( |24l l and, we conclude that 

2XaW-^ Aatn 

2^") “ l + 2XcxW-^a ^ 


l + 2XaW-^a 

2 iA-tAaat£)-i 


X = 2 ( D-t _ 


1 -f 2 Aat£i-ia 
2D~^Xa.cy.^n 


^Dn 


“ l + 2XaW-^a 

Using equation ( |2^ we can express the objective value of 
( fTSl l as 

itr{(x-n)tZl(x-n)} 

tr{ntaatAZl-i(-2)iA(-2)iA-tAaatn} 

“ 2 ( 1 -f 2 Aati:)-ia )2 

tr{ntaatn}4A^(at£)-ia) 

2jY+2XaXD-^ 

2 A^||atn||^(at£)“tQ,) 


( 1 -f 2 Aati:)-ia )2 


= 2 ||uf (at/? la), 


where in the last equality we made use of ( |25] ). We now re¬ 
call that from the third equation in the KKT conditions we 
have that ||u||/A = r -f r' and so we conclude that 


2 (atD-ia) 



jjatnijA 
r + r' J 


and therefore. 


(27) 


itr{(x-n)ti^(x-n)} 

^ ^%(ati:)-ia )2 

(r + r' ± ||a^'n|j)^ 
2 (atZ/”ia) 



{a^D ^a) 
(28) 


Since we are seeking the global minimum, and since we are 
assuming that ||Q:^n|j = \\a{n — n') + {1 — a){fi — n')\\ < 
{r -f r') and hence the constraint is active, we conclude that 


^tr{(x*(a) — n)^ D{x* {a) 

{r + r' — |ja^'n||)^ 1 

2(ati4-ia) ^ 2 


-n)} 


(29) 


Above we have used the fact that = aAn-f (1 — a) An 
and that <yXD~^a = jp -f (1 — ot)^ jp- This proves that 

f is valid when ||a'f'n|j = '^a{n — nl) + {l — a){fL—n')\\ < 
-f r') as long as 1 -f 2Aa'i9“ia = |ja'l'n||/(r -f r') ^ 0. 
Recall that we have already proved that Q holds when when 
||a^'n|| = ||a(n — n') -I- (1 — a)(n — n )|| > {r + r'). 

To finish the proof we now prove the validity of (|^-(|^ 
when ||a^n|| = ||Q!(n — n') + (1 — a){n — n')\\ > {r + 
r'). Recall that we have already proved their validity when 
||a^'n|j = ||a(ri — n!) + {1 — a){n — n')\\ < {r + r'). Now 
notice that, when II a'^njl = ||a(n—n')-|-(l — a)(n — ri')|| > 
(r + r'), we can write. 


2{aW-^a) 

h{a) 



||a'fn|| A 

r + r' J 


2(r + r')\/ cxX D~^(y. 


(30) 


If we define 7 = i_|_ 2 Aat we can write the following 

expression for x*{a) 


x*{a) = n — "fD ^aa^n, (31) 

which holds if ||a^n|j = \\a{n — ri!) + {l — a){n — n')\\ 7 ^ 0 . 
With a little bit of algebra one can see that this is exactly the 
same expression as (|^-(|^ and finish the proof. 





























